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We consider the problem of constructing Wannier functions for Z2 topological insulators in two 
dimensions. It is well known that there is a topological obstruction to the construction of Wannier 
functions for Chern insulators, but it has been unclear whether this is also true for the Z2 case. We 
consider the Kane-Mele tight-binding model, which exhibits both normal (Z2-even) and topological 
(Z2-odd) phases as a function of the model parameters. In the Z2-even phase, the usual projection- 
based scheme can be used to build the Wannier representation. In the Z2-odd phase, we do find a 
topological obstruction, but only if one insists on choosing a gauge that respects the time-reversal 
symmetry, corresponding to Wannier functions that come in time-reversal pairs. If instead we are 
willing to violate this gauge condition, a Wannier representation becomes possible. We present 
an explicit construction of Wannier functions for the Z2-odd phase of the Kane-Mele model via 
a modified projection scheme followed by maximal localization, and confirm that these Wannier 
functions correctly represent the electric polarization and other electronic properties of the insulator. 
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I. INTRODUCTION 

In the past several years there has been a surge of in- 
terest in topological insulators. These are materials that 
are gapped in the bulk, just like ordinary insulators, but 
that cannot be adiabatically connected to ordinary insu- 
lators without closing the gap or breaking some specified 
symmetries. They also exhibit chiral metallic edge states 
that are topologically protected from disorder. 1 3 Topo- 
logical insulators can be distinguished from normal ones 
based on the manner in which the Bloch eigenfunctions 
are topologically twisted in k-space. 

Two types of topological insulators have received the 
most attention. First, Thouless et aL A pointed out long 
ago that a two-dimensional (2D) insulator is character- 
ized in general by a topological integer known as the 
"Chern number" or "TKNN index." A prospective in- 
sulator having a non-zero value of this integer would 
be known as a "Chern" or "quantum anomalous Hall" 
insulator. The latter name arises because such a crys- 
tal would exhibit a quantum Hall effect (QHE) even in 
the absence of a macroscopic magnetic field, and would 
have chiral edge states just like the ordinary field- induced 
QHE. Haldane devised an explicit tight-binding model 
realizing such a case. 5 Since the Hall conductance is odd 
under the time-reversal (T) operator, Chern insulators 
can only be realized in systems with broken T symmetry, 
e.g., insulating ferromagnets. Despite the fact that these 
possibilities have been appreciated now for almost three 
decades, no known experimental realizations of a Chern 
insulator are yet known. 

Second, a great deal of interest has surrounded the 
recent discovery of a different class of topological insu- 
lators known as Z2 insulators that realize the quantum 
spin Hall effect (QSH). 6 Subsequent theoretical 7 9 and 
experimental 10 14 work has succeeded in identifying sev- 
eral materials systems that realize the case of a Z2 topo- 
logical insulator. Unlike the Chern index, which van- 



ishes unless T is broken, the Z2 index (which takes val- 
ues of and 1, or equivalently, "even" and "odd") is 
only well defined when T is conserved. Z2 insulators 
are thus non-magnetic, although a spin-orbit or similar 
interaction is needed to mix the spins in a non-trivial 
way. Because T is preserved, the occupied states at k 
and — k form Kramers pairs, and one can associate a 
Z2 invariant with the way in which these Kramers pairs 
are connected across the Brillouin zone. 15 Since the Z2 
index cannot change along an adiabatic path that is ev- 
erywhere gapped and T-symmetric, a Z2-even (normal) 
insulator cannot be connected to a Z2-odd (topological) 
one by such a path. In 2D there is a single Z2 invari- 
ant, and T-invariant insulators are classified as "even" 
or "odd," while in 3D there are four Z2 invariants and 
the classification is more complicated. 16 

Wannier functions (WFs) have proven to be a valuable 
tool when working with semiconductors and insulators, 
providing a real-space description that can be used to 
understand bonding, construct model Hamiltonians, and 
directly compute certain physical properties such as the 
electric polarization. 17,18 Thus, it is desirable to under- 
stand the construction of the Wannier representation for 
topological insulators so that this useful set of techniques 
can be applied to these novel materials. 

For Chern insulators it has been shown that a non- 
zero Chern number presents a topological obstruction 
that prevents the construction of exponentially localized 
WFs. 19 ' 20 Conversely, a general proof has been given that 
exponentially localized WFs should exist in any 2D or 3D 
insulator having a vanishing Chern index. 21 In principle 
this applies to Z2-odd as well as Z2-even T-invariant insu- 
lators, suggesting that a Wannier representation should 
be possible in both cases. However, it is unclear whether 
the nontrivial topology of the Z2-odd case has any effect 
on the Wannier representation. In particular, one may 
wonder whether the procedure for obtaining WFs would 
be the same as for ordinary insulators, and if not, how it 
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should be modified in order to get well localized WFs in 
the Z2-odd regime. 

In this paper we address this question using the model 
of Kane and Mele 6 as a paradigmatic system that ex- 
hibits both Z2-odd and Z2-even phases. We demonstrate 
that the usual projection scheme used for constructing 
the Wannier representation is still applicable to the Z2- 
odd insulators, but only for gauge choices that do not 
allow WFs to come in time-reversal pairs. We present 
an explicit projection procedure for constructing well- 
localized WFs in the topologically non-trivial phase, and 
show that the WFs can be made even more localized us- 
ing the standard maximal- localization procedure. 17 We 
also discuss the electric polarization from both Berry- 
phase and Wannier points of view, showing the relations 
between the viewpoints and confirming that both give 
identical results. 

The paper is organized as follows. In Sec. II we de- 
fine the Z2 topological invariant in 2D and briefly dis- 
cuss methods for determining it numerically. We review 
the model of Kane and Mele in Sec. Ill, and describe its 
spectrum and phase diagram. In Sec. IV we present the 
projection scheme used to construct WFs and explain 
how the application of this scheme to Z2-odd insulators 
is different than for ordinary insulators. The localiza- 
tion properties of the constructed WFs are described in 
Sec. V. The electric polarization properties and locations 
of the Wannier charge centers are considered in Sec. VI. 
Finally, we make concluding remarks in Sec. VII. 



II. 



Z 2 INVARIANT 



Here we briefly review some of the equivalent ways of 
determining the Z2 invariant in 2D insulators. 

In the work of Ref. 22 the definition of the Z2 invariant 
was given in terms of a function P(k) defined as 



p(k) = pf[( Wi (k)|0K(k))], 



(i) 



i.e., the Pfafflan of a certain k-dependent antisymmetric 
N x TV matrix, where N is the number of occupied bands. 
Here \uj(k)) = e~* k r |'0j(k)) is the periodic part of the 
Bloch function of the j'th occupied band and 6 = is y C 
is the time-reversal operator (C is complex conjugation 
and s y is the second Pauli matrix). If the zeros of P(k) 
are discrete, then the Z2 invariant is odd if the number of 
zeros of the Pfafflan within one half of the Brillouin zone 
(BZ) (see Fig. 1) is odd, and even otherwise. If the zeros 
of the Pfafflan occur along lines in the BZ, then the Z2 
invariant depends similarly on whether half the number 
of sign changes of P(k) along the boundary of the half 
BZ is odd or even. Using A = and 1 to represent 
evenness and oddness respectively, the Z2 invariant can 
equivalently be determined as 6 



A=— / dk • V k log[P(k + iS)] mod 2, 

2lTT J dr 



(2) 



1- 



dr 



r 2 
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FIG. 1. (Color online) Sketch of the Brillouin zone. The 
Berry curvature of Eq. (4) is calculated in the interior of the 
half zone r (dashed region), while the Berry connection is 
evaluated along its boundary dr (arrows indicate direction of 
integration). Time-reversal-invariant points Ti are shown. 



where the loop integral runs along the boundary dr of 
the half BZ, and the 5 term is included for convergence. 

Another approach to the problem of defining A results 
from considerations of "time-reversal polarization." 23 
Here a spin-pumping cycle is considered and it is shown 
that the Z2 index is given by the difference between the 
time-reversal polarizations at the beginning and the mid- 
point of the cycle. This approach leads to the formula 



(-D A =n 



v/detKfOl 



i = i PfKi\)] 



(3) 



where w mn (k) = (u m (—\a)\0\u n (k)) and I\ are the four 
time-reversal invariant points of the BZ (i.e., those for 
which —T{ = + G with G a reciprocal vector). Note 
that the matrix w mn is not the same as that in Eq. (1). 

The definition in Eq. (3) appears to require a knowl- 
edge of the occupied wavefunctions at only four points in 
the BZ, unlike Eq. (2), for which the wavefunctions must 
be known at all points along the boundary of the half BZ. 
However, Eq. (3) is usually not suitable for numerical im- 
plementation in practice, since the sign of the Pfafflan at 
any one of the four points can be flipped by a relabeling 
of the Kramers-degenerate states at that point. To be 
more explicit, there is a "gauge freedom" in the choice of 
states |iz m (k)), corresponding to a k-dependent N x N 
unitary rotation among the occupied states. Eq. (3) is 
only meaningful when a globally smooth gauge choice 
enforces a relation between the labels at the four special 
k-points. 23 This problem may be avoided in the pres- 
ence of some additional symmetry that can be used to 
establish the labels of the bands at these points. For ex- 
ample, in Ref. 9 it is shown how the presence of inversion 
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symmetry allows for a simplified calculation of A from 
Eq. (3). 

In the absence of inversion symmetry, one can use yet 
another definition of the Z2 index taking the form 23 



one can write the lattice definition of the Z2 invariant as 



A = — 

2tt 



AM- 



Udr 



Tdr 



mod 2, 



(4) 



where A = i 



|Vk|^n) is the Berry connection of 



J\f occupied states and T = Vk x A is the corresponding 
Berry curvature. 24 Of course, if A and T are both con- 
structed from a common gauge that is smooth over r, the 
result would vanish by Stokes' theorem. Thus, Eq. (4) 
is only made meaningful by the additional specification 23 
that the boundary integral of A must be calculated using 
a gauge that respects time-reversal symmetry, i.e., 



|^2n-l(-k)) = %2n(k)), 
|^2n(-k)) = -%2n-l(k)>. 



(5) 



For the case of the nontrivial Z2 state, it turns out to be 
impossible to choose a gauge that satisfies both smooth- 
ness over r and the constraint (5) over dr. In other 
words, A=l signals the existence of the topological ob- 
struction. 

To see how this works more explicitly, the contribu- 
tions to the integral of A over dr are illustrated in Fig. 1. 
We choose a gauge that is periodic, \uj(k)) = |uj(k+G)), 
in addition to satisfying Eq. (5). The contributions of the 
top and bottom segments (solid blue arrows in Fig. 1) 
then cancel because they are connected by a recipro- 
cal lattice vector G. Thus, the gauge needs to be fixed 
only along the left and right boundaries (composed of 
red dashed and gray dotted arrows in Fig. 1), which are 
separated by a half reciprocal lattice vector. At each 
of the special points one state from each Kramers- 
degenerate pair is arbitrarily identified as |^2n-i(I\)), 
and the other is constructed via 



\u2 n (r i )) = -o\u 2n - 1 (r i )). 



(6) 



Then we can make an arbitrary gauge choice along the 
remaining portions of the gray dotted arrows in Fig. 1 - 
e.g., accepting the output of some numerical diagonaliza- 
tion procedure. Finally, the gauge should be transferred 
to the dashed-arrow segments using Eq. (5), where k and 
— k belong to the dotted and dashed segments respec- 
tively. 

Eq. (4) can now be evaluated using a uniform dis- 
cretized mesh K covering the region r, with the time- 
reversal constraint applied to the boundary dr as de- 
scribed above. To do so, define the link matrices 
M M)Tim (k) = (ix n (k)|it m (k + s M )) and the unimodular 
link variables L M (k) = det MJ\ det M M |, where k e K 
and Si (S2) is the step of the mesh in the direction 
of the reciprocal lattice vector Gi (G2). By defining 
Ai(k) = logLi(k) and 

F(k) = log[Li(k)L 2 (k + sOL^O* + s 2 )L 2 - 1 (k)], (7) 



1 



]T ^(k)-^F(k) 

.ke<9r kGr 



mod 2. 



(8) 



For a sufficiently fine mesh there will be no ambiguity in 
the branch choice for the complex log in Eq. (7), since 
the argument of the log must approach unity as the mesh 
becomes dense. Moreover, a change in the branch choice 
determining one of the boundary links A s (k) has no effect 
(mod 2) on Eq. (7), since each A s (k) appears twice as a 
result of the gauge- fixing on the boundary. Thus, once 
the mesh is fine enough so that the branch choices in 
Eq. (7) are all unambiguous, Eq. (8) gives A exactly. 25 



III. 



THE KANE-MELE MODEL 



In their remarkable paper introducing a Z2 topologi- 
cal classification to distinguish a QSH (Z2-odd) insula- 
tor from an ordinary (Z2-even) insulator, Kane and Mele 
(KM) 6 also introduced a model tight-binding Hamilto- 
nian that describes a 2D Z2-odd insulator in some of its 
parameter space. In this section we will describe some of 
the properties of the model suggested therein. 

The KM model is a tight-binding model on a honey- 
comb lattice with one spinor orbital per site. The prim- 
itive hexagonal lattice vectors are ai 5 2 = a/2(y/3y ± x) 
and sites A and B are located at = ay/y/3 and 
t# = 2ay/y/3 respectively. The KM Hamiltonian is 

H = tJ2 4 C 3 + * A SO U i3 C i sZc 3 
<ij> <ij> 

+ i\r ^ c\(s x dij) z Cj + X v ^2^ic\ci, (9) 

<ij> i 

where the spin indices have been suppressed on the rais- 
ing and lowering operators, and t is the nearest-neighbor 
hopping amplitude. In the second term, Aso is the 
strength of the spin-orbit interaction acting between sec- 
ond neighbors, with = (2/v // 3)[di x cb] = ±1 depend- 
ing on the relative orientation of the first-neighbor bond 
vectors di and cb encountered by an electron hopping 
from site j to site z, and s z is the z Pauli spin matrix. 
Next, Ar describes the Rashba interaction 26 that couples 
differently oriented first-neighbor spins, with s being the 
vector of Pauli matrices. Finally, X v is the strength of the 
staggered on-site potential, for which & is +1 and —1 on 
A and B sites respectively. Note that the symmetry of 
the problem is lowered significantly compared to an ideal 
honeycomb lattice, since the on-site staggered potential 
makes the A and B sites inequivalent, while the Rashba 
term breaks s z conservation. 

To proceed, we choose the tight-binding basis wave- 
functions to be 



Xj Mr) = (1/V^)$> ik - R ? 

R 



r (r-R-t,), (10) 



di t(l + 2 cos x cos ^/) di2 — 2£cosxsin?/ 

c?2 Aw di5 2Aso(sin2x — 2 sinx cosy) 

d>3 Ar(1 — cos x cosy) c/23 — ARCOsxsiny 

<^4 — Ar sin x sin y c/24 sin x cos y 



TABLE I. (Color online) Nonzero coefficients appearing in 
Eq. (11), using the notation x = k x a/2 and y = y/3k y a/2 (see 
also Fig. 2). 
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FIG. 2. (Color online) Brillouin zone sketched using co- 
ordinates x = k x a/2 and y = V3k y a/2. Primitive re- 
ciprocal lattice vectors Gi = (27r/a)(l, l/s/3) and G2 = 
(2ir/a)(— 1, 1/V$) correspond to gi = (7r, 7r) and g2 = (— 7r, tt) 
respectively. The black rectangle marks the boundary <9£ of 
the zone used for polarization calculations in Sec. VI. 

where a is a spin index, j = {A, B} denotes the atom 
type, tj is a vector that specifies the position of the 
atom in the unit cell, 27 and R is a lattice vector built 
from the primitive lattice vectors ai and &2. This al- 
lows the Hamiltonian to be written as a 4x4 matrix 
Hj a jt a f(k) = (Xjcrkl^lXj'cr'k), which can be cast in terms 
of five Dirac matrices T a and their ten commutators 
T ap = [r« l^]/(2z) as 

5 5 

H(k) = Y / da(k)r a + ^(k)r^ (11) 

a=l a</3=l 

where the Dirac matrices are chosen to be r 1 ' 2 ' 3 ' 4 ' 5 = 
(I®a x ,I®a z ,s x ®a y ,s y ®a y ,s z ®a y ) with the Pauli 
matrices a k and s k acting in sublattice and spin space re- 
spectively. The dependence of the d a and d a p coefficients 
on wavevector is detailed in Table I using the notation 
x = k x a/2 and y = \/3k y a/2, with the relationship of 
these variables to the BZ being sketched in Fig. 2. 

Since, O^O' 1 = T a and OT^Q- 1 = -T af3 , while 
d a (k) = d a (—"k) and d a p(]t) = — d a p(— k), the Hamil- 
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FIG. 3. Phase diagram of the Kane-Mele model for A^/Aso 
> 0. Arrow illustrates a path crossing the phase boundary by 
varying X v while keeping other parameters fixed. 



tonian (9) is time-reversal invariant, i.e., 0H(k)0~ 1 = 
H(—k.). However, it lacks particle- hole symmetry in the 
sense of Refs. (1-3), because of the action of the on-site 
and spin-orbit coupling terms. In the general classifi- 
cation of topological insulators and superconductors, 1 3 
therefore, the Kane-Mele model falls into the All sym- 
plectic symmetry class, which in two dimensions has a Z 2 
classification. This means that by varying parameters of 
the Hamiltonian of Eq. (9) one can switch between Z2- 
odd and Z2-even phases, with the system experiencing a 
gap closure and becoming metallic at the transition from 
one phase to the other. 

For the present purposes we assume Aso > without 
loss of generality. We also fix X v > 0. For this case, 
the transition between Z2-odd and Z2-even phases is ac- 
companied by a gap closure at the K and K' points (the 
zone-boundary points of three- fold symmetry) in the BZ. 
The energy is independent of t at these points, and Aso 
can be used as the energy sca le. The energy gap is the n 
given by |6>/3 - A,/A SO - V( VAso) 2 + 9(A R /A SO ) 2 |, 
leading to the phase diagram shown in Fig. 3. 

Note that when Ar = the model reduces to two in- 
dependent copies of the Haldane model 5 the Z2 invari- 
ant is odd when the Chern numbers are odd, and even 
otherwise. 28 

In what follows we use t as the energy scale and fix 
the values of the other parameters to be Aso A = 0.6 an d 
Ar/£ = 0.5. Varying the third parameter X v /t allows us 
to switch from the Z2-even to the Z2-odd phase. The 
phase transition occurs at \X v \/t ~ 2.93, with the sys- 
tem in the Z2-odd phase for —2.93 < X v /t < 2.93. As 
discussed above, the energy gap closes at the phase tran- 
sition, and remains open in both the Z2-odd and Z2-even 
phases. 
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IV. GAUGE FREEDOM AND WANNIER 
FUNCTIONS 

A. General considerations 

We now consider the problem of constructing Wannier 
functions (WF) for the Kane-Mele model. We emphasize 
that we mean by this a set of localized functions spanning 
the same space as the occupied Bloch bands. Several 
recent papers have discussed the construction of WFs 
for an enlarged subspace including also some unoccupied 
bands for 3D topological insulators such as Bi2Se3, 8,29 in 
which case there is typically no topological obstruction, 
but this is not the context of the present work. 

We start with the general definition of the WF in cell 
R and with band index n in 2D, 



(r|Rn) = W„(r-R) 



A 



(27T) 2 



-ik-R 



^nk(r), 



BZ 



(12) 

where A is the unit cell area and Bloch wavefunctions 
ipnk are assumed to be normalized within the unit cell. 
This definition is not unique; not only is there the usual 
U(l) gauge freedom associated with a k-dependent phase 
twist of each band n, there is more generally a U(M) 
gauge freedom 



l^nk) 



^Z7 mn (k) |^ mk ) 



(13) 



coming from the fact that the M occupied Bloch bands 
can be mixed with each other by a k-dependent U(M) 
transformation. In fact, it is generally necessary to pre- 
mix the Bloch states using this U(N) gauge freedom in or- 
der that the resulting Bloch-like states (and their phases) 
will be smooth functions of k. However, having done so, 
there is still a large gauge freedom associated with the 
application of a subsequent U(Af) gauge rotation that is 
smooth in k. 

This ambiguity in the gauge choice can be removed 
by applying some criterion to the selection of the 
WFs. Since electrons are expected to be localized in 
insulators, 30 a sensible criterion is that of Ref. 17, which 
specifies maximal localization of the WFs in real space. 
In this approach, which we adopt here, one chooses some 
localized trial functions in order to provide a starting 
guess about where the electrons are localized in the unit 
cell, and obtains a fairly well-localized set of WFs by a 
projection procedure to be described shortly. If desired, 
one can follow this with an iterative procedure to make 
the resulting WFs optimally localized. 17 

Consider an insulator with J\f occupied bands. We 
start with a set of M trial states |t;) located in the home 
unit cell, and at each k we project them onto the occupied 
subspace at k to get a set of Bloch-like states 



|T ik ) = P k \ n ] 



N 

E 

n=l 



|^nk)(^nk|Ti). 



(14) 



Since this set of states will not generally be orthonormal, 
we make use of a Lowdin orthonormalization procedure 
which consists of constructing the overlap matrix 

S mn (k) = (T mk |T nk ) (15) 

and obtaining the orthonormal set of Bloch-like orbitals 



/2 



l^mk)- 



(16) 



Note that the -0 n k are not eigenstates of the Hamiltonian, 
but they span the same space, and have the same form, as 
the usual Bloch eigenstates. For an insulator whose gap 
is not too small, and for a set of trial functions embodying 
a reasonable assumption about character of the localized 
electrons, the ^ n k will be smooth functions of k. In that 
case, by the usual properties of Fourier transforms, the 
WFs constructed in analogy with Eq. (12), 



|Rn) 



dkt 



-ik-R 



BZ 



l^nk), 



(17) 



should be well localized. 

Such a construction will break down if the determinant 
of £(k) vanishes at any k. This is guaranteed to occur 
in a Chern insulator, where time-reversal symmetry is 
broken and the Chern index of the occupied manifold is 
non-zero; in this case, construction of exponentially lo- 
calized WFs becomes impossible. 19 21 For a Z 2 insulator, 
however, the presence of time-reversal symmetry guaran- 
tees a zero Chern index, so that exponentially localized 
WFs must exist. 21 In this case, we should be able to find 
a set of trial functions such that det S'(k) ^ throughout 
the BZ. 



B. Z2-even phase 

Let us first apply the method described above to the 
case of the Z2-even phase of the Kane-Mele model. This 
phase is topologically equivalent to the ordinary insula- 
tor, so we anticipate a picture in which the two electrons 
per cell are opposite-spin ones approximately localized 
on the lower-energy (B) site. One way to see this is to 
look at the weights of the basis states in the occupied sub- 
space. Figure 4(a) shows the distribution of these weights 
along a high-symmetry line in the BZ for the Kane-Mele 
model in its Z 2 -even phase. From the figure it is obvious 
that the two basis states on the B site dominate in the 
occupied subspace over the whole BZ. It is then natural 
to choose the two trial functions to be opposite-spin spa- 
tial (5-functions localized on the B site in the home unit 
cell. We choose these to be spin-aligned along z, i.e., 



We?, 



where \af) = |tz) and 
k-space we get 



5(r-t B )|of> (18) 
= \iz)- Transforming to 



Rk) 



R 



e ik - R S(r-K-t B ). 



(19) 



6 



1.0 
0.8 
0.6 
0.4 
0.2 
0.0 

1.0 
0.8 
0.6 
0.4 
0.2 
0.0 



"... 


AT - 
A 1 - 


- bt / ^ M ; 

D v \ p / ~ 






\ / 


r 




K M K" r 
















' "S / 

\ 




-(b) 


\ 

* ♦ \ 

\ 

* - a \ 






FIG. 4. (Color online) Sum of the weights of the projections 
into the two occupied bands of the basis states \A; "fz), \B;^ Z ), 
\A;Xz), and \B;^ Z ) plotted along the diagonal of the BZ for (a) 
X v /t = 5 (Z2-even phase) and (b) X v /t = 1 (Z2-odd phase). 
Inset in (a): BZ of a honeycomb lattice. 

The two occupied Bloch bands may be written as 

l^nk) =Y, C ^\Xi\.) (20) 

i 

where t is a combined index for sublattice and spin, t = 
{1, 2, 3, 4} = {A t, B t, A |, 5 I}, and x^k = Xjak are the 
tight-binding basis functions of Eq. (10). With Eq. (19) 
the projected functions become 

|Tik> = C 2 * lk |Vik> + C 2 * 2k |V>2k>, (21) 

|T 2k )=C 4 * lk |^ lk )+C 4 * 2k |V2k). (22) 

The overlap matrix S is constructed from these functions, 
and for the determinant one finds 

det [S(k)} = (|C 21k | 2 + |C 22k | 2 )(|C 4 i k | 2 + |C 42 k| 2 ) 

— l^lkC^lk + ^22kC'42k| 2 - (23) 

Recall that for the Lowdin orthonormalization procedure 
to succeed, this determinant must remain non-zero ev- 
erywhere in the BZ. This is indeed the case for the Z 2 - 
even phase, as illustrated in Fig. 5(a), where the solid 
black curve shows the dependence of the determinant on 
k along the high-symmetry line in the BZ. 

In contrast, the dashed red curve in Fig. 5(a) shows 
the behavior of det [*S(k)] in the Z 2 -odd regime. The de- 
terminant can be seen to vanish at the K and K' points 
in the BZ. Clearly, this choice of trial functions is not ap- 
propriate for building the Wannier representation in the 
Z 2 -odd phase. Indeed, as we shall see in the next subsec- 
tion, any choice of trial functions that come in Kramers 
pairs is guaranteed to fail in the Z 2 -odd case. There we 
shall also investigate alternative choices of trial functions 
that allow for a successful construction of WFs. 



FIG. 5. (Color online) Plot of det[5(k)] along the diagonal of 
the BZ for X v /t = 5 (Z2-even phase) and X v /t = 1 (Z2-odd 
phase), (a) Trial functions are \B;^ Z ) and \B;^ Z ). (b) Trial 
functions are |A;tx) and \B\^ X ). 

C. Z2-odd phase 

To gain some insight into the appropriate choice of trial 
functions in the Z 2 -odd regime, consider the weights of 
the basis functions in the occupied space shown for this 
case in Fig. 4(b). Unlike the normal insulator, the Z 2 -odd 
phase does not favor any particular basis states. Instead, 
different basis states dominate in different portions of the 
BZ. For example, at points K and K' the occupied space 
is represented by only two of the four basis states; at each 
of these points the two participating basis states have 
opposite spin and sublattice indices, and none appear 
in common at both points. (The states at K are, of 
course, Kramers pairs of those at K' .) It follows that 
if any of the trial states is simply set equal to one of 
the four basis states, then at least one of the |T) would 
vanish either at K or K' , and the determinant would 
vanish there too. This explains the failure of the naive 
Wannier construction procedure for the Z 2 -odd phase; 
with the naive choice of trial functions as in Eq. (18), the 
determinant vanishes at both K and K', as shown by the 
red dashed curve in Fig. 5(a). 31 

In fact, this failure can be understood from a general 
point of view. If the two trial functions form a Kramers 
pair, then the projection procedure of Eqs. (14-16) will 
result in Bloch-like functions obeying 

|^i(-k)) = 6# 2 (k)), 
hM-k)) = -%ii(k)). (24) 

The WFs obtained from Eq. (17) will then also form a 
Kramers pair. But Eq. (24) is nothing other than the 
constraint of Eq. (5) denning a gauge that respects time- 
reversal symmetry, and it has been shown 23,32 ' 33 that an 
odd value of the Z 2 invariant presents an obstruction 
against constructing such a gauge. In other words, in 



7 



the Z 2 -odd phase a smooth gauge cannot be fixed by 
choosing trial functions that are time-reversal pairs of 
each other, and a choice of WFs as time-reversal pairs is 
not possible. Hence, in order to construct the Wannier 
representation in the Z 2 -odd regime, one should choose 
trial functions that do not transform into one another 
under time reversal. 

Following these arguments, we choose the two trial 
functions to be localized on different sites in the home 
unit cell. Moreover, in order that they will have compo- 
nents on states with spins both up and down along z, we 
choose the spins of the trial states so that one is along 
+x and the other along — x. M In k-space this becomes 

i r *) = ^=E elk ' R<5 ( r - R -^) ( 25 ) 



where ti = and t 2 = t^, leading to 



|Tik> = [(C? lk + C 3 * lk )|Vi> + (Cr 2k + C* 32k )\ip 2 )} /y/2 

(26) 

and 



|T 2k ) = [(C 2 * lk - C 4 * lk )|Vi> + (C 2 * 2k - C 4 * 2k )|^ 2 )] /y/2. 

(27) 

The determinant takes the form 



det[5] = (|Cii k + C 3 i k | 2 + |C 12k + C 32k | 2 )(|C 21k - C 4 i k | 2 + |C 22k - C 42k | 2 )/4 - 

- |(Cnk + C f 3ik)(C2i k - CJik) + (Ci2k + ^32k)(C 22k - C4 2k )| 2 /4. (28) 



The dependence det[*S(k)] is shown along the diagonal 
of the Brillouin zone for this choice of trial functions in 
Fig. 5(b). In the Z 2 -odd phase (dashed line) the determi- 
nant remains non-zero everywhere in the BZ. 35 Not sur- 
prisingly, the same trial functions are very poorly suited 
to the normal-insulator phase, as can be seen from solid 
line in the same panel. In this case det[S r (k)] almost van- 
ishes at K and K' and remains quite small throughout 
the rest of the BZ, so that one should clearly revert to 
the time-reversed pair of trial functions of Eq. (18) and 
Fig. 5(a) in order to get well-localized WFs. 



We made an arbitrary choice above in selecting the two 
trial functions to be up and down along x. In fact, if we 
repeat the entire procedure using trial functions that are 
spin-up and spin-down along any unit vector ft lying in 
the #2/-plane, we find that det[5(k)] changes very little, 
with only small changes in the size of the dip near the Y 
point. Thus, we find that the choice of trial functions in 
Eq. (25) is not unique. Instead, there is a large degree of 
arbitrariness in the choice of WFs in the Z 2 -odd case. 



To conclude, we have established that the choice of a 
time-reversal pair of trial functions, Eq. (18), that al- 
lows for the construction of well-localized WFs in the 
ordinary-insulator phase cannot be used in the Z 2 -odd 
phase. In order for the usual projection method for 
constructing the Wannier representation to work in this 
topologically nontrivial phase, the trial functions should 
explicitly break time-reversal symmetry, i.e., they should 
not come in time-reversal pairs. 



V. LOCALIZATION OF WANNIER 
FUNCTIONS IN THE Z 2 -ODD INSULATOR 

Now that we know how to construct WFs for the Z 2 - 
odd insulator, we discuss their localization properties. As 
we have noted in the preceding section, the choice of the 
trial functions, Eq. (25), is not unique; there are other 
gauge choices arising from different trial functions that 
also produce well-defined sets of WFs. Since different 
gauge choices lead to different degrees of localization of 
the resulting WFs, it is natural to fix the gauge by the 
condition of maximal possible localization of the WFs. 

The problem of constructing maximally-localized WFs 
was studied by Marzari and Vanderbilt. 17 They consid- 
ered the total quadratic spread 



O = X)[<0n|' 



2 1 On) - (On|r|On) 2 



(29) 



as a measure of the derealization of WFs in real 
space, and developed methods for iteratively reducing the 
spread via a series of unitary transformations, Eq. (13), 
applied prior to WF construction. The spread functional 
was decomposed into two parts, Q = Qj + O, with 



E 

n=l 



(On\r 2 \On) - ^ |(Rra|r|On)| 2 



Ftm 



being the gauge- invariant part and 



^ = E E m™\r\On)f 

n=l Rm^On 



(30) 



(31) 



the gauge-dependent part of the spread. Discretized k- 
space formulas for Eqs. (30) and (31) were also derived 
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FIG. 6. (Color online) Wannier spreads Qi and Cl for the 
Kane-Mele model on a 60x60 k-mesh, initialized using the 
trial functions of Eq. (25). "Initial" and "final" values are 
those computed before and after the iterative minimization 
respectively. The system is in the Z2-odd phase for X v /t < 
2.93. 

for the case that the BZ is represented by a uniform k 
mesh. The resulting expression for the gauge- invariant 
spread is, for example, 

n ' = jf E w * £ - i^ k+b) i 2 ) . (32) 

k,b m,n=l 

where 

4 

^i k „ k+b) = («nk|n mk+b ) = ^c; nk Q mk+be - ibt ^ 

<?=1 

(33) 

are overlap matrices and b are "mesh vectors" connect- 
ing each k-point to its nearest neighbors. The latter are 
chosen, together with a set of weights cj^, in such a way 
as to satisfy the condition 

^2uj b bibj = 5ij. (34) 

b 

A corresponding expression for and a description of 
steepest-descent methods capable of minimizing ^, were 
also given in Ref. 17. Note that, in order to avoid getting 
trapped in false local minima, the iterative procedure 
is normally initialized using the trial-function projection 
procedure described in Sec. IV above. 

We now apply this method to the Kane-Mele model. 
The lattice is hexagonal, and in this case six bj vectors 
are needed to satisfy the condition (34), namely bi = 

bi = Gi/q, b 2 = b~, (Gi + G 2 )/q, and b 3 
— b6 = G 2 /q. All six have the same length b and weight 
= 1 / (36 2 ). We start with the WFs obtained with the 
projection method using the trial functions of Eq. (25), 
appropriate for the Z2-odd phase. 

The resulting spreads, both before and after the it- 
erative minimization, are shown in Fig. 6. (f2j, being 
gauge- invariant, is the same before and after.) The left 
part of the figure shows the behavior in the Z2-odd phase, 
where the trail functions are the appropriate ones. The 
results in this region were not strongly sensitive to the 



k-point mesh density. The fact that Q is similar in mag- 
nitude to the unminimized Qi, and that the localization 
procedure reduces ^ by only 20 — 30%, provide additional 
evidence that the choice of trial functions was a good one. 
The Wannier charge centers were almost unchanged by 
the minimization procedure; the ^-coordinates were zero, 
while fiy ~ a/y/3 and f 2y — 2a/ a/3 (see Sec. VI for 
details), in good agreement with our initial assumption 
about the WFs being localized on A and B sites. 

The right part of Fig. 6, for X v /t > 2.93, shows what 
happens when we attempt to use the same trial functions 
in the normal phase. Qi is of course unaffected by the 
choice of trial functions, and the fact that it has a smaller 
value in this region indicates, not surprisingly, that the 
insulating state is simpler and more localized in the nor- 
mal state. (For large X v /t the WFs approach spatial 
delta functions, explaining the fact that asymptotes 
to zero in that limit.) Not surprisingly, however, using 
the trial functions appropriate to the Z2-odd phase in the 
Z2-even regime results in very poor localization of the 
WFs as measured by Q. Our data also suggests that in 
the Z 2 -odd phase MLWFs are less localized than MLWS 
in the Z2-even phase. For example, the use of trial func- 
tions (18) with X v /t = 5 and a 60 x 60 k-mesh results in 
ft/ = 0.02770 and Q = 0.00025. We also find that the 
results are more sensitive to the choice of k-mesh in the 
Z2-odd regime. 

To summarize the results of this section, we studied the 
construction of maximally localized WFs in the Z2-odd 
phase using the Kane-Mele model as an example. We 
have seen that our initial guess of Sec. IV about the lo- 
calization of WFs in this topological regime is very good, 
and that the maximal localization procedure does not 
greatly reduce the spread. 

VI. HYBRID WANNIER CHARGE CENTERS 
AND POLARIZATION 

In this section we discuss the polarization in Z 2 -odd 
insulators using the example of the Kane-Mele model, 
and see what insights about the topological insulating 
phase can be obtained by inspecting this property. 

The electronic polarization in a 2D system can be de- 
fined either in terms of the Berry phase 36 

II N r 

P = #Wlm^ / dk(*i nk |V k K k ) (35) 

( Z7V) n=l J 

or via the summation of Wannier charge centers 18 

P = -^i>' ( 36 ) 

where e is the electronic charge and A is the area of the 
unit cell. The two definitions are identical and define 
electronic polarization modulo a polarization quantum 
|e|R/A, R being a lattice vector. This ambiguity can 
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be understood as a freedom in the choice of branch in 
Eq. (35) or in the choice of unit cell in Eq. (36). The 
definition via Wannier charge centers makes the depen- 
dence of P on the choice of origin obvious. As described 
in Sec. Ill, the origin of the Kane-Mele model is chosen 
such that atoms are located along the y-ax.is at = £y/3 
and t# = 2£y/3, where £ = |ai + a2| = ay/3. Because the 
Hamiltonian has 3-fold symmetry, we expect the rescaled 
polarization (A/|e|)P to lie at the origin, at t^, or at t^- 
To distinguish between these possibilities it is sufficient 
to compute P y , which is well-defined modulo \e\/a. 



A. Total polarization 

A direct computation of electronic polarization via 
Eq. (35) in the Z2-even phase results in P y = |e|/3a mod 
|e|/a, consistent with the fact that both Wannier cen- 
ters in Eq. (36) lie at t B (since -4\e\£/3A = — 8|e|/3a = 
|e|/3a mod \e\/a.) In the Z2-odd phase, on the other 
hand, Eqs. (35) and (36) lead to P y = mod \e\/a. 
Again, this is consistent with the locations of the WFs. 
As indicated in Sec. V, the Wannier centers f n in this 
phase lie approximately at and t#. More precisely, 
we find that they are located at f i = (1 — 5)£y/3 and 
?2 = (2 + S)£y/3, where S is a small correction (e.g., 
S = 0.0018 at X v /t = 1). Thus, the sum of the Wannier 
centers is just £y, or zero modulo a lattice vector. 

It is interesting to note that, in retrospect, the compu- 
tation of the polarization via Eq. (35) would have given 
a strong hint about the appropriate choice of trial func- 
tions in the Z2-odd insulator. That is, knowing only that 
P y = 0, one might have guessed that both WFs should be 
centered halfway between and t#, or both at the cen- 
ter of the honeycomb ring, or one at and the other at 
t#. The latter possibility becomes the most likely when 
we also take into account that in the Z2-odd phase the 
two WFs cannot form a Kramers pair. 



B. Hybrid Wannier decomposition 

In order to obtain a deeper understanding of the ori- 
gin of the polarization and expose some qualitative differ- 
ences in the behavior of its k-dependent decomposition in 
Z2-even and odd phases, it is useful to use a hybrid repre- 
sentation in which the Wannier transformation is carried 
out in one direction only. As indicated above, we know 
from symmetry considerations that we can set P x = 
and characterize the polarization by P y mod £|e|/A To 
compute P y , it is convenient to choose the BZ to be a 
rectangle extending over k x G [0, 2tt / 'a) and k y G [0, 47r/£] 
(corresponding to the region ( in Fig. 2). We can then 
define hybrid WFs 

t />47r/<e 

\nkj y ) = ^J dky |^ nk ) (37) 



in terms of which the usual WFs are 

|Rn) = \nl x l y ) = ^- [ ' dk x e- ik * l *\nk x l y ). (38) 

The hybrid Wannier centers are defined as 

Vn(k x ) = (nk x 0\y\nk x 0) (39) 

and the total electronic polarization is 

IpI _ r 27r / a 
Py = -— C Y J dk x y n (k x ). (40) 

In practice the k x integral is discretized by a sum over a 
mesh of k x values, and at each k x the y n (k x ) are calcu- 
lated by considering the corresponding string of k-points 
along k y . In the case that the gauge has been specified 
by a particular set of 2D WFs |Rn), or, equivalently, 
by the corresponding Bloch-like functions l^nk), this is 
done straightforwardly using the discretized Berry-phase 
formula 

VM = -^lm\og]jM^ (41) 

3 

where is a shorthand for the overlap matrix 

M^^j+i) of Eq. (33) connecting /^-points j and j + 1 
along the string. 

As was emphasized in Sec. IV, the ^ n k carry the in- 
formation about the gauge choice. Thus, different gauge 
choices - i.e., different choices of WFs - will result in 
different hybrid WFs and different y n (k x ). However, the 
sum y n (k x ) at a given k x is gauge- invariant, and as a 
result P y of Eq. (40) must remain the same in any gauge. 

Of special interest is a gauge choice in which, at each 
k x , the hybrid WFs \nk x l y ) are maximally localized in 
the y direction. It was shown in Ref. 17 that in ID the 
Wannier charge centers can be obtained by a parallel- 
transport construction using the overlap matrices M^ J \ 
Specifically, the "unitary part" of each overlap ma- 

trix is obtained by carrying out the singular-value de- 
composition M = VT>W\ where V and W are unitary 
and E is real-positive and diagonal, and then setting 
M = VW^ . This is reasonable because, for a sufficiently 
fine mesh spacing, E is almost the unit matrix. Then, the 
unitary matrix A = . describes the transport of 
states along the string. The eigenvalues A n of this matrix 
are all of unit modulus, and their phases define Wannier 
centers via 37 

Vn (k x ) = - — Im log A n . (42) 

47T 

Note that no iterative procedure is needed. Inserting this 
equation into Eq. (40), one gets a discretized formula for 
P y that is consistent with Eq. (35). 



10 





3 










2 




- (a) 






1 






n 







k x a 


2n 




k x a 


2n 


















2 








: (b) 


s, - - — -*^/^=~=-- 




1 


:(d) 






: ' ■ ■—^7 


n 





FIG. 7. (Color online) Hybrid Wannier centers y n (k x ), in 
units of £/2, for the Kane-Mele model. Z2-even phase (X v /t — 
3): (a) maxloc gauge; (b) WF gauge of Eq. (18). Z2-odd phase 
(X v /t — 1): (c) maxloc gauge; (d) WF gauge of Eq. (25). In 
each case, several periodic images are shown. 



C. Results 

We illustrate these ideas now for the KM model in its 
normal and Z2-odd phases. In each case we present re- 
sults for y n (k x ) for two choices of gauge: the maximally- 
localized one along y as discussed in the previous para- 
graph, and the one corresponding to the WFs constructed 
from the trial functions of Eq. (18) for the Z2-even phase 
or those of Eq. (25) for the Z2-odd phase. In what fol- 
lows, we refer to these as the "maxloc" and "WF-based" 
gauges respectively. 

In the ordinary insulating regime, the maxloc and 
WF-based y n (k x ) curves look very similar to each other. 
Fig. 7(a) and (b) show the calculated results for the case 
of X v /t = 3, very close to the transition on the insu- 
lating side (recall the critical value is at X v /t = 2.93). 
Three of the infinite number of periodic images along y 
are shown. The "bumps" in the curves near the K and 
K' points in the BZ are the result of the proximity to the 
transition; as one goes deeper into the insulating phase, 
the curves flatten out and become smooth functions of 
k x . The solid and dashed curves are mirror images of 
each other; in the maxloc construction of Fig. 7(a) this 
just reflects the time-reversal invariance of the Hamilto- 
nian, while in Fig. 7(b) it follows from the fact that the 
WFs form a Kramers pair. 

When averaged over k x , each curve is found to have 
a mean y value of 2£/3 to numerical precision, or £/6 
modulo £/2, consistent with the discussion in Sec. VIA. 

The corresponding results for the Z2-odd phase are 
shown in Fig. 7(c) and (d) for X v /t = 1. As expected, 
there is again a mirror symmetry visible in the curves 
for the maxloc construction in Fig. 7(c), but the con- 
nectivity of the curves is qualitatively different: in going 
from k x = to n/a we see that the n'th solid curve 
goes up to cross the (n + l)'th dashed curve, while the 
n'th dashed curve goes down to cross the (n — l)'th solid 
curve. This is exactly the kind of behavior that was ex- 
hibited in Fig. 3(a) of Ref. 23 as a signal of the Z 2 -odd 



phase. Moreover, if we follow the n'th dashed curve all 
the way across the BZ, we find that it wraps to become 
the (n + l)'th one when k x = 2tt / 'a wraps back to k x = 0. 
This is precisely the kind of behavior that is characteris- 
tic of a Chern (or quantum anomalous Hall) insulator, 38 
which implies that we can assign a Chern number of +1 
to the Bloch subspace spanned by the eigenvectors cor- 
responding to the dashed bands. However, since we are 
studying here a system with time-reversal symmetry, we 
find also a partner subspace corresponding to the full 
curve in Fig. 7(c) having Chern number —1. As a result, 
of course, the overall occupied space has a total vanishing 
Chern number, as it must due to the time-reversal sym- 
metry. The evaluation of the polarization P y through 
Eq. (40) again yields P y = mod \e\/a, consistent with 
the direct calculation of Sec. VIA. 



Finally, Fig. 7(d) shows the y n {k x ) curves for the same 
Z2-odd parameters as in Fig. 7(c), but using the WF- 
based gauge determined by the trial functions of Eq. (25). 
At any given k x , we confirm that y\ + yi is the same in 
Fig. 7(d) as in Fig. 7(c), and the total polarization is 
therefore the same. However, because the two WFs do 
not form a Kramers pair in this case, the dashed and solid 
curves do not map into each other under time-reversal 
symmetry, and there is no degeneracy at k x = ir/a. 
Moreover, the Chern number of each band is individu- 
ally zero, consistent with the fact that each one is de- 
rived from a WF. The average y values for the solid and 
dashed curves are 0.978^/3 and 2.022^/3 mod f/2, very 
close to the nominal locations of the trial functions at 
and t#, respectively. 



To recap, in both the Z2-even and Z2-odd cases, we 
find that the occupied Bloch space can be cast as the 
direct sum of two subspaces that map into one another 
under the time-reversal operation, corresponding to the 
solid and dashed curves of Figs. 7(a-c). These subspaces 
are not built from Hamiltonian eigenstates, but from suit- 
able k-dependent U(2) rotations among the Hamiltonian 
eigenstates. In the Z2-even case the Chern index of each 
of these subspaces is separately zero, so that we can also 
provide a Wannier representation for each subspace sepa- 
rately. This is essentially the case of Fig. 7(b), and since 
the spaces form a time-reversal pair, the WFs form a 
time-reversal pair as well. In contrast, for the Z2-odd 
phase, the decomposition into two subspaces that are 
time-reversal images of each other necessarily results in 
subspaces having individual Chern numbers of ±1, and 
these are not individually Wannier-representable. Only 
by violating the condition that the two spaces be time- 
reversal partners, as was done in Fig. 7(d), can we de- 
compose the space into two subspaces having zero Chern 
indices individually. By doing so, we can find a Wannier 
representation of the entire space, but only on condition 
that the two WFs do not form a Kramers pair. 
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VII. CONCLUSIONS 



In this paper we have considered the question of how 
to construct a Wannier representation for Z2-odd topo- 
logical insulators in 2D. We have shown that the usual 
method based on projection onto trial functions fails be- 
cause of a topological obstruction if one imposes the 
condition that the trial functions should come in time- 
reversal pairs. On the other hand, the projection method 
can be made to work if this condition is not imposed, re- 
sulting in WFs that do not transform into one another 
under time reversal. 

Such a Wannier representation may have some formal 
disadvantages. For example, if one writes the Hamil- 
tonian as a matrix in this Wannier representation, its 
time-reversal invariance is no longer transparent, and the 
presence of other symmetries may become less obvious 
as well. On the other hand, it does satisfy all the usual 
properties of a Wannier representation, as for example 
the ability to express the electric polarization in terms of 
the locations of the Wannier centers, and there is every 
reason to expect that the maximally localized WFs are 
still exponentially localized. 21 

The generalization of our findings to the 3D case should 
be relatively straightforward. Certainly the topological 
obstruction to the construction of Kramers-pair WFs re- 



mains for both weak and strong Z 2 topological insulators 
in 3D. To see this, consider in turn each of the six symme- 
try planes in k-space (fci = 0, k<i = 0, ks =0, k\ = tt/cl, 
etc.) on which behaves like a 2D time-reversal in- 
variant system. For both weak and strong topological 
insulators, at least one of these six planes must have a 
Z2-odd 2D invariant. But if a gauge exists obeying the 
time-reversal condition of Eq. (5) in the 3D k-space, then 
it does so in particular on the 2D plane, in contradiction 
with the 2D arguments about a topological construction. 

Thus, the general strategy for constructing WFs for 3D 
topological insulators should be very similar to the one 
presented here in 2D. Namely, one has to choose pairs 
of trial functions that do not transform into one another 
by time-reversal symmetry, and to do it in such a way 
that the projection of these trial functions onto the Bloch 
states does not become singular anywhere in the 3D BZ. 
While it may be interesting to explore how this might 
best be done in practice for real 3D topological insulators, 
e.g., in the density- functional context, the choice is likely 
to depend sensitively on details of the particular system 
of interest. Thus, an investigation of these issues falls 
beyond the scope of the present work. 
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